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Abstract 

A class of smoothing methods is proposed for solving mathematical programs with 
equimibrium constraints. We introduce new and very simple regularizations of the 
complementarity constraints. Some estimate distance to optimal solution and expan- 
sions of the optimal value function are presented. Limited numerical experiments 
using SNOPT algorithm are presented to verify the efficiency of our approach. 

1 Introduction 

Mathematical programs with equilibrium constraints (MPECs) constitute an important 
class of optimization problems and pose special theoretical and numerical challenges. 
MPECs are constrained optimization problems in which the essential constraints are de- 
fined by some parametric variational inequalities or a parametric complementarity system. 
MPECs can be closely related to the well-known Stackelberg game and to general bilevel 
programming. As a result, MPECs play a very important role in many fields such as en- 
gineering design, economic equilibrium, multilevel game, and mathematical programming 
theory itself, and it has been receiving much attention in the optimization world. 
However, MPECs are very difficult to deal with because, the feasible region and optimal 
solution set are almost non convex non concave and not even connected. Moreover, the 
constraints can not satisfy any standard constraint qualification such as the linear inde- 
pendence constraint qualification or the Mangasarian-Fromovitz constraint qualification 
at any feasible point [U [13] . 

In this paper, we consider MPECs in their standard complementarity constrained opti- 
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mization problems formulation 

min f(x, y) 

s.t. x e X,y g R m ,z e R l , A G R l , 

F(x,y)-V y g(x,y) T X = (1.1) 
c/(x,y) = z 
z > 0,A > 0,A T z = 

where the functions / : R n+m -> E, F : R n+m -> M m and 5 : M n+m -> IR' are all twice 
continuously differential) le and X is a nonempty and compact subset of M". 



Remark. The constraints of ( 1.1 ) correspond to the KKT conditions of the parametrized 
variational inequality 

y G C{x) and (v - y) T F(x,y) > for all v G C(x), (1.2) 

where C(x) := {y G R m /5(x,y) > 0}. 

The negative properties of MPECs make these problems very difficult and exclude any 
direct use of standard non linear programming (NLP) algorithms. 

In this paper we propose some smoothing techniques to regularize the complementarity 
constraints and construct relaxed problems that are suitable for NLP algorithms. 
Many regularization and relaxation techniques have already been proposed, here is an 
incomplete list of such methods 
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£ln{e t 


+ 6^1 = Vi 



(1.3) 

In almost all these techniques, the constraints XiZ{ = or min(Aj, zi) = are replaced by 
some smooth approximations. 

In our approach, the complementarity constraint is interpreted componnent-wise as: 

Vi, At most one of z% or Aj is nonzero. 
So, we construct some parameterized real functions that satisfy: 



(0 r (a?) ~ 1 if x / 0) and (9 r (x) ~ if x = 0) 
to count nonzeros and then replace the constraint 
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by 

9 r (\i) + e r ( Zi ) < i. 



In section 2, we present some preliminaries and assumptions on the problem (1.1) (es- 
sentially the same as in [TT]) . In Section 3, the smoothing functions and techniques 
are presented and many approximation and regularity properties are proved. Section 4 
is devoted to the analysis of the regularization process. The last section presents some 
numerical experiments concerning two smoothing functions. 



2 Assumptions and preliminaries 



We essentially need the same assumptions and background as in [llj. A complete presen- 
tation of this background needs about 6 to 7 pages. We will only present in this section 
some definitions, known optimality conditions and constraint qualifications. For some 
others we will only refer readers to |llj . These notions will be useful in the next section. 
The first definition concern a first order optimality condition: the strong stationarity 



Definition 2.1 A feasible point (x* , y* , z* , X*) is strongly stationary for (1.1) if d = 

solves 



min Vf\x*,y*Y d XiV 

s.t. d x G Add(X{x*)),d z eR l ,d x eR l , 

VF(x*,y*) T d x , y - V y g(x*,y*) T d x - V(V y g(x*, y*)) T d x , y = 



d z = 



(2.1) 



Vg(x*,y*yd x 
{d z )i =0,i G I z \h 
(d x )i = 0, i G I\\I Z 
{d z h>o,(d x )i>o,ie i z ni x 

where d = (d x , d y , d z , d x ) T G R n + m + 2i 5 I z and I x are the active sets at (x* ,y* , z* , X*) 

I g :={i = l,..., l\z* = 0} and I x := {t = 1, . . . , l\ X* = 0} 

and Add{X{x*)) is the admissible directions set defined by 

Add{X{x*)) := {d x G R m \3r > V0 < r < r x* + rd x G X .} 

Remark. There is an other kind of stationarity (the B-stationarity) which is less restric- 
tive but very difficult to chek. We prefer to not present it in this paper. These two 
stationarity properties are equivalent when the MPEC-LICQ (defined next) is satisfied 
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Definition 2.2 The MPEC-LICQ is satisfied at the point (x* , y*, z*, A*) if the linear in- 
dependence constraint qualification (LICQ) is satisfied for the following RNLP problem at 
(x*,y*,z*,X*). 

min f(x, y) 

s.t. x£X,z£R l ,\£R l , 

F{x,y)-X7 y g{x,y) T \ = 

g(x,y) = z (2.2) 
Zi = 0,i £ I Z \I\ 
Xi = 0,i£ h\Iz 

zi > o, \i > o,i e i z n i x 

V 

An other important and usefull constraint qualification is the following Mangasarian- 
Fromovitz one 

Definition 2.3 The MPEC-MFCQ is satisfied at the point (x* , y* , z* , A*) if the Mangasarian- 
Fromovitz constraint qualification (MFCQ) is satisfied for the RNLP problem at (x*,y*, z* , A*). 

We will also use some Second-Order sufficient conditions namely: the (MPEC-SOSC) 
and the (RNLP-SOSC). These two conditions (among others) are defined in 

3 The smoothing technique 

For r > 0, we consider real functions 9 r : K+ — ► [0, 1] satisfying 

(i) r is nondecr easing, strictly concave and continuously differentiable, 
(it) Vr > 0, r (O) = 0, 

(3.1) 

(Hi) Vx > 0, ]im9 r (x) = 1, and 

r— +0 

(iv) lim ^(0) > 0. 

We will present some interesting examples of such functions after the following approxi- 
mation result 

Lemma 3.1 For any e > 0, and x,y > 0, there exists ro > such that 

Vr < ro, (min(x,y) = 0) =>■ (0 r (x) + 9 r (y) < 1) =>■ (min(x,y) < e). 

Proof - The first property is obvious since 9 r (0) = and r < 1. 
Using assumption (Hi) for x = e, we have 

Va > 0, 3r > 0/ Vr < r 1 - 9 r (e) < a, 
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so that, if we suppose that min(x,y) > e, assumption (i) gives 

r (x) + r (y) > 20 r (e) > 2(1 - a). 

Then if we choose a < h, we obtain that 9 r (x) + 9 r (y) > 1. □ 
This first approximation result can be improved for some interesting choices of the smooth- 
ing functions 9 r 

m (*) = ^ 

{9 Wk ) 0W»(x) = l-e-W foik>0 

nlogx „loe/_\ log(l+x) 



log(l + x + r) 

We will also consider the general class O- 1 of functions 



5>1 



verifying (i — iv) and 



3>1 



> 



Remark. The functions 9 Wk are the density functions of Weibull distributions, when 
k = 1, the obtained smoothing method corresponds (with slight modifications) to the 
inequality entropic regularization [2]. Simple comparison calculus proove that #! og and 
9 Wk for (0 < k < 1) belong to the class of functions O- 1 . 



Lemma 3.2 we have 

(i) Vx > 0, Vy > 

(ii) Vx > 0, Vy > 

Proof - (i) We have 



0*0*0 + 0r(y) ^ 1 < ^ = ^ > x ■ y < r 2 , and 
x • y = = 



' r >l(s) + 0£l(y) <1 



x ■ y <r l 



l 1 r {x) + 9 l r {y) 



so that 



r 1 (x) + r 1 (y) < 1 



2xy + rx + ry 
xy + rx + ry + r 2 ' 

2xy + rx + ry < xy + rx + ry + r 2 
x • y < r 2 . 



The first part of (ii) follows obviously from Lemma 3.1 and the second one is a direct 
consequence of (i) since 

0}\x) + 9}\y) < 1 =► 9}.{x) + 9 l r {y) < 1. 

□ 



Using any function 9 r satisfying (|3.1|), we obtain the relaxed following problem for (1.1 ) 



(3.2) 



min f(x,y) 

s.t. (x,y) G G E^, A G R l + ,e G 
F(x,y)- Vj / y(x,y) T A = 
y(x,y) = z 

0r(Ai) + r (zi) + f% = 1, Vi G {1, I}. 
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Remarks, (i) By choosing some particular smoothing functions (ex. 0^ k ), the nonnega- 
tivity constraints on A and z become implicite and can be removed from the definition of 
( |3.2[ ).(This can have an important impact in practice.) 

(ii) Under some classical assumptions, as in [5] we can easily prove that the jacobian of 
equality constraints (with respect to (y, z, A) is nonsingular. This property is useful in 
practice since standard NLP algorithms use Newton-type to solve systems of nonlinear 
equations corresponding to this jacobian. 

Lemma 3.3 // g is concave with respect to y and F is uniformly strongly monotone with 
respect to y, then for every nonnegative r and every feasible point (x, y, z, A, e) of problem 



(3.2), the jacobian of equality constraints (with respect to (y,z,X)) is nonsingular. 



Proof - Using the assumptions 3.1 (i) and (iv), the proof is exactely the same as in [2] or 
0. □ 



Problem (3.2) may be viewed as a perturbation of (1.1). Previous lemmas prove that 



(3.2) is in fact some tight relaxation of (1.1). However this perturbation is not continuous 
on the parameter r so that any direct use of perturbation results such that [3] is impossi- 
ble. 

Fortunately, Lemma 3.2 proves that for the particular smoothing function 9}, the corre- 
sponding relaxed problem (3.2) is equivalent to (Reg(r 2 )) in [11]. We can then benefit 
from the theoretical results in [llj . 

The following results provide, in the case of the 6} function, some distance estimate be- 



tween solution of (3.2) and solution of (1.1). These results correspond to applications of 
[[3], Theorem 5.57, Theorem 4.55 and Lemma 4.57] and can be found with complete proofs 
in [TT]. We just state them in our context and add the optimal value expansion. 



Theorem 3.1 Suppose that X* = (x* ,y* , z* , A*) is a strongly stationary point of ( 1 . 1 ) at 
which MPEC-MFCQ and MPEC-SOSC are satisfied. Then there are positive constants 
a, f, and M such that for all r S (0,f], the global solution X{r) of the localized problem 
\3.2fy with the additional ball constraint \\X — X*\\ < a that lies closest to X* satisfies 



\\X(r) — X*\\ < M.r. Furthermore the optimal value v r of (3.2) has an expansion of the 
form 

v r 

where v° is the optimal value of (1.1) and a is the optimal value of an auxiliary quadratic 
problem^. 



n 2 / 2 \ 

v + — .a.r + o(r J 



Theorem 3.2 Suppose that X* = (x*, y* , z* , A*) is a strongly stationary point of (1.1) 
at which MPEC-LICQ and RNLP-SOSC are satisfied. Then there are positive constants 
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a, r, and M such that for all r 6 (0,r], the global solution X{r) of the localized problem 
\3.2fy with the additional ball constraint \\X — X*\\ < a that lies closest to X* satisfies 



\\X(r) — X*\\ < M.r 2 . Furthermore the optimal value v r of (3.2) has an expansion of the 
form 

v r < v° + b.r 2 + 0(r 4 ) 



where v is the optimal value of (1.1) and b is the optimal value of an auxiliary linearized 
problem^. 



For functions of the general class 0- , the corresponding feasible sets satisfy 



J~P C T ' n>\ C J~Q 



where Tp, F e >i and T e \ are respectively the feasible set of problem (1.1) and (3.2) for the 
corresponding 9 r function. 



These inclusions prove that the optimal value expansions given in Theorem 3.1 and The- 
orem 



3.2 are still valid under the same assumptions. 



Theorem 3.3 When using functions Or 1 , under the same assumptions of Theoren 3.1 



(resp. Theoren 3.2) the optimal value v r of (3.2) has an expansion of the form 



I 

v r < v + - .a.r 2 + o(r 2 ) (resp. v r < v° + b.r 2 + 0(r 4 )) 



4 Numerical resuts 

For two different smoothing functions, we present some numerical results using the SNOPT 
[8] nonlinear programming algorithm on the AMPL [1] optimization plateform. Our aim 
is just to verify the qualitative numerical efficiency of our approach. We consider a subset 
of the MACMPEC [9] test problems with known optimal values and solutions (a large 
part of these test problems were used by [2j [5] in their numerical experiments) . 
We choose the two functions 



X 



x + r 



and 



9^(x) = l-e 



The first function has (in our analysis) the best theoretical results and corresponds "in 
some way "to the regularization studied in |12[ [TT] . While the second one corresponds to 
the enropic regularization [2J |6] . 
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In our experiments, we made a logarithmic scaling for these two functions to bound their 
gradients. Each constraint 

r (Ai) + Or(zi) + e { = l 
is in fact replaced by the following inequality 

r 2 ln ( + > 0, 

\\i + r Zi + rJ 

in the case of the Q\ function and 

e » • + e I > 0. 

in the case of the function. 

The two following tables give for each considered problem and for different starting points, 
the used value of the parameter r, the optimal value and solution obtained when using 
each of the two smoothing functions. The tables report also different informations con- 
cerning the computational effort of the solver SNOPT. itM and itm correspond to the 
total number of major and minor iterations numbers |8j. The total number of objective 
function evaluations is given in (Obj.). {grad.) corresponds to the total number of objec- 
tive function gradient evaluations, (constr.) and (jac.) give respectively the total number 
of constraints and constraints gradient evaluations. 
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Problem 


r 


Start 


Obj.val. 


Opt.x 


(itM.itm) 


Obj. 


grad 


constr. 


Jac 


Bardl 


l.e-2 


no 


17 


(1,0) 


(5,8) 


9 


8 


9 


8 


Dfl 


l.e-3 


no 





(1,0) 


(1,1) 


3 


2 


3 


2 


Gauvin 


l.e-2 


no 


20 


(2,14) 


(4,H) 


7 


6 


7 


6 


jrl 


l.e-2 


no 


0.5 


(0.5,0.5) 


(6,3) 


9 


8 


9 


8 


GnashlO 


l.e-5 


gnashl0.dat 


-230.8232 


47.036 


(17,46) 


21 


20 


21 


20 


Gnashll 


l.e-4 


gnashll.dat 


-129.9119 


34.9942 


(20,50) 


18 


17 


21 


20 


Gnashl2 


l.c-4 


gnashl2.dat 


-36.93311 


18.1332 


(24,51) 


27 


26 


27 


26 


Gnashl3 


l.e-2 


gnashl3.dat 


-7.061783 


7.55197 


(14,56) 


20 


19 


23 


22 


Gnashl4 


l.e-3 


gnashl4.dat 


-0.179046 


1.06632 


(14,46) 


18 


17 


21 


20 


Scholtesl 


l.e-1 


1 


2 





(9,10) 


14 


13 


14 


13 


Bilevell 


l.e-2 


(25,25) 


5 


(25,30) 


(3,11) 








9 


8 






(50,50) 


5 


(25,30) 


(0,6) 








2 


1 


Nashl 


l.e-1 


(0,0) 


1.61e-14 


(9.996,4.999) 


(13,42) 


25 


24 


25 


24 






(5,5) 


1.60e-18 


(9.313,5.686) 


(10,33) 


32 


31 


32 


31 






(10,10) 


1.46e-14 


(9.092,5.901) 


(16,38) 


34 


33 


34 


33 






(10,0) 


3.56e-24 


(9.999,4.999) 


(12,34) 


28 


27 


28 


27 






(0,10) 


9.03e-22 


(9.999,4.999) 


(14,41) 


31 


30 


31 


30 


Bilevel2 


l.e-4 


(0,0,0,0) 


-6600 


(6.441,4.863,12.559,16.137) 


(6,43) 


9 


8 


9 


8 






(0,5,0,20) 


-6600 


(6.575,5,12.425,16) 


(6,50) 


9 


8 


9 


8 






(5,0,15,10) 


-6600 


(6.837,12.162,16) 


(5,36) 


7 


6 


7 


6 






(5,5,15,15) 


-6600 


( 4.892,3.373,14.107,17.627) 


(5,35) 


7 


6 


7 


6 






(10,5,15,10) 


-6600 


(8.014,4.971,10.986,16.029) 


(5,38) 


7 


6 


7 


6 


Bilcvcl3 


l.e-4 


(0,0) 


-12.6787 


(0,2) 


(9,23) 


12 


11 


12 


11 






(0,2) 


-12.6787 


(0,2) 


(15,27) 


32 


31 


32 


31 






(2,0) 


-10.36 


(2,0) 


(01,06) 


3 


2 


3 


2 


dcsilva 


l.e-3 


(0,0) 


-1 


(0.5,0.5) 


(4,11) 


6 


5 


6 


5 






(2,2) 


-1 


(0.5,0.5) 


(3,9) 


5 


4 


5 


4 


Stack. 1 


l.e-2 





-3266.6666 


93.3333 


(4,9) 


6 


5 


6 


5 






100 


-3266.6666 


93.3333 


(3,3) 


5 


4 


5 


4 






200 


-3266.6666 


93.3333 


(7,5) 


19 


18 


19 


18 



Tablel: using the 9} smoothing function 
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Problem 


r 


Start 


Obj.val. 


Opt.x 


(itM,itm) 


Obj. 


grad 


constr. 


Jac 


Bardl 


l.e-2 


no 


17 


(1,0) 


(13,8) 


16 


15 


16 


15 


Dfl 


l.e-3 


no 





(1,0) 


(1,2) 


3 


2 


3 


2 


Gauvin 


l.e-2 


no 


20 


(2,14) 


(5,12) 


7 


6 


7 


6 


jrl 


l.e-2 


no 


0.5 


(0.5,0.5) 


(13,4) 


16 


15 


16 


15 


GnashlO 


l.e-3 


gnashl0.dat 


-230.8232 


47.036 


(17,63) 


19 


18 


21 


20 


Gnashll 


l.e-3 


gnashll.dat 


-129.9119 


34.9942 


(15,48) 


18 


17 


18 


17 


Gnashl2 


l.c-1 


gnashl2.dat 


-36.93311 


18.1332 


(15,43) 


19 


18 


19 


18 


Gnashl3 


l.e-1 


gnashl3.dat 


-7.061783 


7.55197 


(23,69) 


30 


29 


30 


29 


Gnashl4 


l.e-3 


gnashl4.dat 


-0.179046 


1.06633 


(22,38) 


27 


26 


27 


26 


Scholtosl 


l.e-1 


1 


2 





(11,11) 


16 


15 


16 


15 


Bilevell 


l.e-2 


(25,25) 


5 


(25,30) 


(3,11) 








9 


8 






(50,50) 


5 


(25,30) 


(0,6) 








2 


1 


Nashl 


l.e-1 


(0,0) 


7.27e-14 


(9,6) 


(9,16) 


12 


11 


12 


11 






(5,5) 


4.25e-18 


(10,5) 


(6,16) 


10 


9 


10 


9 






(10,10) 


1.09e-ll 


(9,6) 


(13,25) 


21 


20 


21 


20 






(10,0) 


1.27e-13 


(9.355,5.645) 


(16,34) 


27 


26 


27 


26 






(0,10) 


3.29e-15 


(9.396,5.604) 


(7,16) 


11 


10 


11 


10 


Bilevel2 


l.e-1 


(0,0,0,0) 


-6600 


( 4.851,5,14.149,16) 


(6,34) 


8 


7 


8 


7 






(0,5,0,20) 


-6600 


(5.195,5,13.805,16) 


(5,40) 


7 


6 


7 


6 






(5,0,15,10) 


-6600 


(6.099,4.834,12.901,16.166) 


(5,42) 


7 


6 


7 


6 






(5,5,15,15) 


-6600 


(4,1.714,15,19.286) 


(5,45) 


7 


6 


7 


6 






(10,5,15,10) 


-6600 


(7.724,5,11.276,16) 


(5,50) 


7 


6 


7 


6 


Bilevel3 


l.e-1 


(0,0) 


-12.6787 


(0,2) 


(22,38) 


29 


28 


29 


28 






(0,2) 


-12.6787 


(0,2) 


(34,61) 


56 


55 


56 


55 






(2=0) 


-10.36 


(2,0) 


(01,06) 


3 


2 


3 


2 


dcsilva 


l.e-2 


(0,0) 


-1 


(0.5,0.5) 


(5,9) 


8 


7 


8 


7 






(2,2) 


-1 


(0.5,0.5) 


(6,14) 


9 


8 


9 


8 


Stack. 1 


l.e-2 





-3266.6666 


93.3333 


(4,4) 


6 


5 


6 


5 






100 


-3266.6666 


93.3333 


(3,5) 


5 


4 


5 


4 






200 


-3266.6666 


93.3333 


(11,4) 


14 


13 


14 


13 



Table2: using the inequality enropic approach (O^ 1 ) 



5 Conclusion 

We introduced a new regularization scheme for mathematical programs with complemen- 
tarity constrains. Our approach is very simple and quite different from existing techniques 
for the same class of problems. The obtained regularized problems are now suitable for 
standard NLP algorithms. These regularizations have different theoretical sensivity and 
regularity properties. The limited numerical experiments give very promising results (com- 
parable to those of [2]) and suggest to make real investigations on functions of the class 
9 Wk . Therefore, we hope that some of our smoothing functions will correspond to simple 
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and efficient algorithms for the solution of real- world MPECs and Bilevel programs. 
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